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ABSTRACT 

Cosmic rays are a fundamental source of ionization for molecular and diffuse clouds, 
influencing their chemical, thermal, and dynamical evolution. The amount of cosmic 
rays inside a cloud also determines the y-ray flux produced by hadronic collisions be¬ 
tween cosmic rays and cloud material. We study the spectrum of cosmic rays inside 
and outside of a diffuse cloud, by solving the stationary transport equation for cos¬ 
mic rays including diffusion, advection and energy losses due to ionization of neutral 
hydrogen atoms. We found that the cosmic ray spectrum inside a diffuse cloud differs 
from the one in the interstellar medium (ISM) for energies smaller than E^r ~ 100 
MeV, irrespective of the model details. Below the spectrum is harder (softer) 
than that in the ISM if the latter is a power law oc with s larger (smaller) than 
- 0.42. 
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1 INTRODUCTION 


The amount of penetration of cosmic rays (CRs) into molec¬ 
ular clouds (MCs) regulates the ionization leve l of clouds 
and dense cores (fo r reviews see lDaIgarn3 l2006l : ICeccarellil 
I2OIII : Ilndriololl2013l ) and thus affects their dynamical evolu¬ 
tion and the process of star formation. Moreover, the exclu- 
sion of CRs from MCs can reduce their gamm a-ray emission 
dSkilling &: StrongIl976l : fGabici et al.ll2007t) . which results 
from the decay of neutral pions produ ced in inelas tic in¬ 
teractions of CRs in the dense gas (see lGabicill201 j . for a 
review). Therefore, it is of prime importance to understand 
wether CRs do penetrate or not MCs. 

The difficulty of modeling the problem of CR penetra¬ 
tion into clouds resides in its highly non-linear nature: CRs 
generate magnetic turbulence at the cloud border due to 
streaming instability. The level of the turbulence in turn 
determines the diffusion coefficient of CRs and thus, pre¬ 
sumably, their capability of penetrate clouds. Some dis¬ 
crepancy exists in the literature, different theoretical ap¬ 
proaches to the problem giving different results. According 
to early papers, CRs with energies below t ens or hundreds of 
MeV are effectively excluded from MCs dSkilling &: Strong 
19761: ICesarskv fc Voli3ll978l ). while in a more recent work. 


Everett fc Zweibell ( 20 111) found out that the CR intensity is 


only slightly reduced inside clouds. However, a direct com¬ 
parison of the two approaches is not straightforward, since 
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the former are kinetic approaches, while the latter a two 
fluid ones. An implicit assumption in all these papers is the 
fact that streaming instability would enhance the magnetic 
turbulence and cause the exclusion of CRs from clouds. 

Here, we present a solution of the steady-state ki¬ 
netic transport equation of CRs along a magnetic flux tube 
that encompasses a MC. We generalize the simplified two- 
zo nes (in- and out-side o f the cloud) kinetic approac hes 
bv ISkilling fc Strong lll976l ) and lCesarskv fc V5U3 lll978l ) by 
considering the full spatially dependent equation. Remark¬ 
ably, we find that the exclusion of CRs from diffuse clouds 
of typical column density Nh « 3 x 10^^ cm“^ is effective 
below an energy of ~ 100 MeV, independently on the pres¬ 
ence or not of streaming instability. In fact, the exclusion 
energy Ef,r depends only (and quite weakly) on the phys¬ 
ical parameters that characterize the ISM and the gas in 
the cloud. This result suggests that: i) the suppression of 
the gamma-ray emission from a cloud due to CR exclusion 
is not significant (the threshold for neutral pion production 
being ~ 280 MeV), and ii) the intensity of CRs is suppressed 
inside MCs at the particle energies which are most relevant 
for ionization. 


2 APPROXIMATE ANALYTIC SOLUTION 

Consider a cloud of size Lc and of hydrogen density uh 
threaded in a magnetic field of intensity Bq, oriented along 
the x-axis. Such a one-dimensional configuration is realis¬ 
tic if one considers spatial scales smaller than the mag- 
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Figure 1. Sketch of the simplified 1-D model used to describe 
the cloud geometry. 


netic field coherence length in the ISM, i.e. ~ 50-100 pc. 
Moreover, the magnetic field strength is assumed to be 
spatially constant and not to change across the transition 
from outside to inside of the cloud. This is supported by 
observations (ICrutcheil 120101 1 showing that the magnetic 
held strength is independent on the density of the inter¬ 
stellar medium as long as the latter remains smaller than 
« 300 cm“®. The cloud is assumed to be immersed in a 
diff use, hot, and fully ioniz e d ISM of density Uj. F ollow¬ 
ing ISkill ingfc St rong! dlOTOl l , ICesarskv fc Volkl dlOTSl l , and 
lEverett fc Zweibell (l201ll l. CRs are assumed to propagate 
along the magnetic held lines only, i.e. diffusion perpendic¬ 
ular to held lines is set to zero. 

For simplicity, the transition between the low density 
and ionized ISM and the dense cloud is taken to be sharp 
and located ai x = Xc- In order to determine the CR prohle 
at such transition, three regions are dehned (see Figure [T]): 
(1) a zone far away from the cloud, dehned as a: < 0, where 
the CR density is unaffected by the presence of the cloud, 
is spatially homogeneous and equal to the sea of Galactic 
CRs, described by the particle distribution function fo{p) oc 
p”®, (2) a zone immediately outside of the cloud (0 < a: < 
Xc) where the CR density f{p, x) is expected to be space 
dependent and to differ from /o(p) due to the presence of 
the cloud itself, and (3) the dense and neutral intra cloud 
medium, dehned as Xc < x < Xc + Lc, with a spatially 
averaged CR density of {fc{p))- 

Due to severe ionization energy losses in the dense 
cloud, one expects {fc{p)) < /o(p), and thus a negative 
spatial gradient of the CR density forms in region (2). A 
gradient in the CR density drives the streaming instability, 
characterized by a growth rate of Alfven waves equal to : 

_ 16 tt^vaVpP'^ df _ rgjR df 

3B^W{kr) dx~ W dx ’ ^ ’ 


where va is the Alfven speed, kr ~ l/rg(p) the resonant 
wavenumber for particles of Larmor radius Vg, and Vp and p 
the particle velocity and momentum, respectively. The spec¬ 
tral distribution of Alfven waves is described by W{k), de¬ 
hned as (SB)^ = B^ f W{k)dk/k, and is believed to regulate 
the penetration of CRs into the cloud by determining the 
value of the CR diffusion coefficient: 


37vW{kr) 


( 2 ) 


The streaming instability excites Alfven waves that 
move with velocity va towards the MC. Thus, if D 2 is 


the value of the diffusion coefficient in zone (2), assumed 
to be constant in space, dimensional analysis suggests that 
Xc ~ D 2 Iva- This means that CRs that try to escape diffu¬ 
sively from the MC are advected back to it by Alfven waves 
if they are within a distance Xc from the MC border. The 
steady-state CR transport equation in zone (2) then reads: 

VAd^f = D2dlf (3) 


where we have neglected the energy loss term which is unim¬ 
portant in low density environments. This equation can be 
solved in an approximate way by reminding that region (2) is 
characterized by the condition x < D 2 /UA, i.e. the diffusion 
term dominates over the advection one. If we then drop the 
advection term VAdf /dx the solution of Equation [3] reads: 

f =-[fo - fix/:)]x/xc +fo (4) 

where f{Xc) is the CR distribution function immediately 
outside of the MC. This approximate solution can now be 
used to compute the flux of CRs into the cloud as: 

- + 2f{xZ)vA = 2foVA (5) 

where the factor of 2 takes into account the fact that CRs 
penetrate the cloud from two sides. This can be rephrased 
by stating that the flux of CRs entering the cloud is inde¬ 
pendent on the value of the diffusion coefficient, and also on 
any detail of the magnetic field amplification (we never used 
Equations [ 1 ] and [ 2 ] to derive it!). The only assumption made 
is that of the existence of Alfven waves converging towards 
the MC. Indeed, such a situation is expected also in the ab¬ 
sence of streaming instability, because waves are damped in 
the dense and neutral gas of the MC, and thus in the cloud 
vicinity one does not expect any appreciable flux of waves 
coming from the cloud. Thus, Equation [5] is valid indepen¬ 
dently on the effectiveness of the streaming instability, and 
the results presented here are very general. 

This is a quite remarkable result, and implies a sort 
of universality of the solution of the problem. Indeed, the 
spectrum of CRs inside the cloud {fc(p)) can be obtained by 
balancing the flux of CRs entering the cloud with the rate at 
which CRs are removed from the cloud due to energy losses 
dSkilling fc Stron'3ll976ll : 

2}o{p)va = [ pp ^(/ c ( p )}] • ( 6 ) 

Equation |6] is valid only when ionization losses play a role (if 
they can be neglected one gets the trivial solution {fc{p)) ~ 
fo{p)), and dimensional analysis suggests that this indeed 
happens when: 


v{p) 


VATljp) 

Lc/2 


(7) 


where ti = —p/p is the characteristic momentum loss time 
due to ionization. This fact can be interpreted as follows: 
since inside the MC waves are strongly damped, CRs are 
expected to free stream at a velocity ~ Vp, where Vp is the 
velocity corresponding to a momentum p. The cloud crossing 
time is thus Tc ~ Lc/vp, and Equation [7] implies that CRs, 
in order to be significantly affected by energy losses, have to 
repeatedly cross the MC, a number of times of the order of 
xi/tc ^ Vp/vA, which can be very large. As said above, for 
T] > 1 the trivial solution is found, {fc{p)) = fo{p), which 
means that the rate at which CRs of a given momentum are 
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removed from the cloud due to ionization losses « foLc/n 
is smaller than the rate at which CRs are advected into the 
cloud ~ 2foVA, i.e. rj > 1. On the other hand, when rj < 1 
the advective flux of CRs into the MC is not sufficient to 
balance the loss rate of CRs due to energy losses. 

In the energy range between 100 keV and 1 GeV the 
loss time n due to ionization of neutral hydrogen can be well 
approximated by a power law in momentum, which reads: 



where the normalization and the slope are tq = 1.46 • 10® yr 
and a = 2.58, respectively, an d have been obt a ined f itting 
the energy losses provided by IPadovani et al.l (120091 ') (see 
their Figure 7). By using this expression we can convert 
Equation [7] into an energy range: 

/ s -0.78 / jy \ O-'^S 

^ - ^%T00^j ( 3xlO^fcm-0 

(9) 

where Nh is a typical column density for a diffuse cloud. 
The expression for the break energy Ebr highlights the main 
result of this paper: the exclusion of CRs from MCs does 
not depend, as one might have expected, on the level of 
magnetic turbulence at the cloud boundary (the diffusion 
coefficient D 2 ), but depends only on the properties of the 
ISM (namely, the Alfven speed in the diffuse phase and the 
column density of the dense phase). This is because the flux 
of CRs into the MC is fixed to « 2 /oua by the self regulating 
spatial gradient of CRs that forms in region (2). 

Above the energy Ebr (or momentum pbr) {fc{p)) = 
fo{p), while below that energy an approximate solution can 
be obtained after integrating Equation [6] 


(Mp)) 


foiPbr) 



X 


ji _ (jlS] 

( S 3 y Pbr J 




which implies that, for p ^ pbr the solution is a power law 
{fc{p)) oc for s < a—3 and {fc{p)) oc for s > 0—3. 
Remarkably, for a CR spectrum in the ISM /o(p) oc p~’‘ with 
s = Q —3 ~ 0.42 the slope of the spectrum of CR is identical 
inside and outside of the cloud. 


3 FORMAL SOLUTION 


In this Section we provide a formal solution of the one di¬ 
mensional steady state equation for the transport of CRs in 
presence of diffusion, advection and energy losses. Written 
in the rest frame of the plasma the equation reads: 

dx 

[?5(a:,p)//(a^,P)] =0, (11) 


D{x,p) 


df{x,p) 


dx 


ox 3 ax Op 


where f{x,p) is the CR distribution function, D{x,p) is 
the diffusion coefficient, p{x,p) is the momentum losses and 
va{x) = B{xf l^/A-Kpi is the Alfven speed in the regular 
magnetic field B(x). Notice that in this expression we in¬ 
clude only the ions mass density, pi, because at the wave¬ 
lengths relevant in this work, the wave frequency is smaller 


than the charge-exchange frequency between ions and neu¬ 
trals, hence, while ions oscillate with waves, neutrals h ave no 
time to couple with them ('e.g. IZweibel fc ShulllIlQ^ i. Note 
that Equation imposes a diffusive behavior to CRs both 
outside and inside of the cloud. The validity of such approx¬ 
imation inside the cloud is questionable, because the mag¬ 
netic turbulence is damped very efficiently by the ion-neutral 
friction. Nevertheless for the moment we assume that the 
propagation is diffusive also inside the cloud. We will show 
that our prediction on the CR spectrum is not strongly af¬ 
fected by this assumption 

Let us now discuss the boundary conditions. Far 
from the cloud {x Xc) we impose that the CR 

distribution reduces to the Galactic one fo{p), while 
at the cloud centre we impose the symmetry condition 
dxf(x,p)\r.=r.. -i-r ,„/2 ~ 0. This is differ ent from what has 
been done in lEverett fc Zweibdl (1201 1[) .where a condition 
on the CR gradient far from the cloud, rather than at the 
cloud center was imposed, to match the CR gradient that 
one would expect in a Galactic environment. However, such 
a choice may lead to an unphysical solution, because it is 
what happens inside the cloud that determines the value of 
the gradient outside of the cloud, and not viceversa. For this 
reason we decided to use the symmetry condition, which is 
valid as long as the diffusion approximation holds. In our ap¬ 
proach the value of dxf away from the cloud is an outcome 
of the calculation. We also notice that the symmetry condi¬ 
tion is valid only for isolated clouds and breaks if the cloud is 
located near a CR source. Such a situation requires different 
boundary conditions and will be considered elsewhere. 

If we introduce the function g — Ddxf, Equation 1111 
reduces to a linear differential equation of the first order: 

dxQ - gvA/D + Q = 0 , (12) 


The nonlinearity of the problem has been hidden in the func¬ 
tion Q{x,p) which plays the role of a source/sink term: 


Q{x,p) 


p dvA df 1 d r . 2 r-i 
3 dx dp dp 


(13) 


The solution for g is: 


nXc-\-Lc/' 2 . 

g{x,p)= Q{x',p)exp 

J X 

(14) 

We can now write the solution above for the simpli¬ 
fied geometry of the MC sketched in Figure [T] Inside the 
MC we assume a constant density of neutral hydrogen of 
riH = 100 cm“®, with a ionization fraction of 10~^. For the 
diffuse ISM we take rn = 10“^ cm~®, so that the Alfven 
speed is constant across the transition between the diffuse 
medium and the MC. For the moment we ignore the ef¬ 
fect of streaming instability and we assume a Kolmogorov 
diffusion coefficient outside of the MC: D{x,p) — Dkoi{p) ~ 
10^^ (p/mcf^^P cm^/s, with j3 = Vp/c. Inside the MC Alfven 
waves are heavily damped due to ion-neutral friction and the 
CR diffusion coefficient is S> D^oi- 

Under these assumptions it is straightforward to derive 
from Equation M an expression for / outside of the cloud: 


f 


VA 


D{y,p) 


dy 


dx'. 


f{x,p) = /o(p)-l-—e' 
VA 


L 


x:c + Lc/2 v^(x'-xc) 

Q{x ,p)e dx . 

( 15 ) 
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which tells that the CR density outside of the cloud is af¬ 
fected up to a distance of: 


Dkoi 

VA 


300 P 


(b^)’ 




(16) 

which can be much larger than the cloud size, and, for a 
particle energy of ~ 100 MeV is of the order of 100 pc, 
comparable to the magnetic field’s coherence length in the 
interstellar medium. For distances to the MC larger than 
Xc the CR distribution reduces to the Galactic one. Strictly 
speaking Equation 1161 implies that, for E > 100 MeV the 
1-D approximation breaks down and a more complex trans¬ 
port model should be adopted. However, when the effect of 
the streaming instability is taken into account (see the 
validity of the 1-D approach is guaranteed up to particle en¬ 
ergies of few hundreds of MeV, well above the critical energy 
defined in Equation [9] 

Equation [15] can be further simplified by assuming that 
the diffusion coefficient inside the cloud is Dc S> LcVa (a 
condition which is easily fulfilled). This implies that a critical 
momentum exists given by the condition that the loss time 
is longer than the CR propagation time across the cloud 
Ti{p*) > Lc/2Dc{p*). Then, for p ^ p* , the CR spatial 
distribution inside the MC is roughly constant and after 
some manipulations Equation 1151 becomes: 

clLc 1 9 


f{x,p) = fo{p) - —e' [pp^fc 

VA 2 p^ dp 


(17) 


where we also made use of the continuity of the CR distri¬ 
bution function at the MC border = f{x'^) = /(a;“). 

We note that Equation [T7] can be used to find the solu¬ 
tion of the problem also when free streaming of CRs (instead 
of diffusion) is assumed inside the MC. This is because also 
under such assumption, a momentum p* exists above which 
the CR distribution function can be considered spatially 
constant. In this case p* is determined by Ti{p*) = Lcjvat 
where Vst ~ Vp/3 is the free streaming velocity of CRs in¬ 
side the MC. For Lc = 10 pc and uh = 100 cm“® we get 
E{p*) ~ 3 MeV. The only difference with respect to the dif¬ 
fusive case is that under the assumption of free streaming 
inside the MC the continuity condition f{xZ) = f{xt) has 
to be considered only as an approximate one. 

After using Equation [8] and performing the derivative 
in momentum. Equation El becomes: 


f{x,p) = fo{p)+ vip)[{^-a)fc+pdpfc]e'' (18) 


and can be solved for x = Xc to give: 

Mp)= [ »7(p)/o(p)exp|- f {a - S + V{P))^\Y 

J p K J p P ) P 

(19) 

As stressed in Section 2, a remarkable property of this solu¬ 
tion is that /c does not depend on the diffusion coefficient 
but only on the Alfven speed and on the MC properties. On 
the other hand, it is easy to show that the CR spectrum 
outside the MC does depend on the diffusion coefficient as: 

fix,p) = fo{p) + [/c(p) - /o(p)] . (20) 

which is shown in Figure [2] together with the approximate 
solution given in EquationjT] for a spectrum of Galactic GRs 
/o(p) oc p~^'’^ normalized to an energy density of 1 eV cm“®. 
Figure [ 3 ] shows few examples for /c(p) obtained from 



Figure 2. Spatial profile of the CR density outside of the MC 
according to the formal (solid) and approximate (dashed) solu¬ 
tion. The dot-dashed line accounts for the presence of streaming 
instability. 



Figure 3. Spectra of CRs inside the cloud (solid lines) assuming 
that the spectra far for the cloud are given by power law in mo¬ 
mentum of slope s (dashed lines). The dot-dashed line shows the 
eigenfunction of the problem. 


Equation [ 19 ] using a simple power law for /o oc p~'‘ (showed 
with dashed lines). We notice that the function /o oc = 
p-o i '2 jg eigenfunction of the problem, giving fc = fo- 


4 STREAMING INSTABILITY 

The main result of this paper is the fact that the CR spec¬ 
trum inside of the MC does not depend on the diffusion co¬ 
efficient in the vicinity of the cloud. However, streaming in¬ 
stability may affect the solution outside of the cloud, by de¬ 
termining the spectral and spatial distribution of CRs there 
fEquations ID and 1201 explicitly depend on the CR diffusion 
coefficient outside of the MC). The full non-linear problem 
is described by two co upled equations and one f or CRs, one 
for Alfven waves (e.g. iLaeage fc Cesarskvlll983l ) 

VAda:f = da:{DB/Wda:f) ( 21 ) 

VAda^W = FcrW - FaW (22) 
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Figure 4. Spatial distribution of the energy density of Alfven 
waves. The lines refer to waves resonating with particles of energy 
E as indicated by the labels. 

where Fd is a damping rate for Alfven waves and Db = 
ivrgl^Tv is the Bohm diffusion coefficient. In the following, 
we will distinguish among the CR diffusion coefficient in the 
ISM Dkoi = Db /Wo and the amplified one in the vicinity of 
the MC, Db/W , where Wo is the typical spectrum of Alfven 
waves in the ISM. 

To estimate the maximum possible effect of streaming 
instability one can set the damping rate to zero Fd = 0. 
Though this is a quite unrealistic as sumption (many pro- 
cesses of wave damping exist, see e.g. iFelice fcKulsrudll200ll 
and references therein), it nevertheless allows to estimate a 
strict lower limit for the CR diffusion coefficient in the vicin¬ 
ity of a MC. If Fd = 0, Eauations l21landl22lhave an analyti c 
solution (see Equations 8 and 9 in lLagage fc Cesarskvlfl98^ : 

f = fo- 7-r-^- (23) 

(1 - ) exp [-(a: - xPj/xc\ - 1 

Wo 

W = Wo + - -r-^- (24) 

- jf^) exp [-(* - xP)/xc\ - 1 

where /c = fix~), a = F)s/(FcrXc) and Xc = Duoi/va- 
For X = x'^ the waves grow up to Wc = W(Xc) = 
Wo + (rcR/uA)[/o — f{xp)]. The extension of the region Ax 
affected by the presence of the cloud can be estimated by 
imposing in Equation 1231 f = (/o -I- fc)/2 which gives Ax = 
a;cln[(Wo + IFc)/II4], which tends to zero for Wc ^ Wo and 
to Xcln(2) « Xc for Wc = Wo- For example, if streaming in¬ 
stability increases the energy density of the turbulent field 
by an order of magnitude, then Ax ~ Q.lXc- This implies 
that for significant field amplification the extension of the 
region affected by the presence of the MC is strongly reduced 
with respect to the estimate given in Eguation llbl which was 
derived from the assumption of no amplification. Thus, in 
this case the 1-D approximation holds, since we deal with 
distances much shorter than the field coherence length. 

A more quantitative estimate is provided in Figures [2] 
and m where Equations [23] and (231 have been plotted by 
assuming for /o(p) the interstellar spectrum of CRs (s = 4.7 
and a total energy density of 1 eV cm“®). The dot-dashed 
curve in Figure[3]refers to Alfven waves resonating with CRs 
of particle energy equal to Ebr, while the solid one to waves 
resonating with particles of energy 1 MeV. An amplification 


of the turbulent field energy density W of roughly an order 
of magnitude is achieved over a region of size <C Xc- 

To conclude, it is easy to check that also in the presence 
of streaming instability the flux of CRs entering the MC is 
2foVA- 


5 DISCUSSION AND CONCLUSIONS 

In this paper we demonstrated that the flux of CRs enter¬ 
ing a MC is of the order of 2foVA- This result holds both 
in the presence or absence of turbulent magnetic field am¬ 
plification due to CR streaming instability. Due to the bal¬ 
ance between advective flux of CRs into the MC and energy 
losses into the cloud, an equilibrium spectrum forms inside 
the MC, characterized by a feature at an energy of the order 
of Ebr ~ 100 MeV for diffuse MCs. Below Ebr the spectrum 
falls below (rises above) the CR spectrum in the ISM if the 
latter is steeper (harder) than oc This fact will have 

an impac t on the estimates of the CR ionization rates in 
MCs fsee lPadovani et al.ll2009l . for a review). However, be¬ 
fore drawing firm conclusions on this issue, the role of CR 
electrons must be assessed, and this will be done in a forth¬ 
coming publication. 

The results obtained here are based on the assumption 
of stationarity. This assumption is justified because the typ¬ 
ical time scale of the problem can be estimated from dimen¬ 
sional analysis as D/va which for a particle energy of 100 
MeV gives ~ 10^{D/Dkoi){v a/WO km/s)“^ yr, which is al¬ 
ways shorter than the dynamical (free-fall) time of the cloud 
Rs 10’^(n_R/100 cm“®) yr. 

Finally, we note that a break in the spectrum in the 
100 MeV range would not affect significantly the gamma- 
ray luminosity of the cloud, because the threshold for neutral 
pion production in proton-proton interactions is at a larger 
energy, namely « 280 MeV. 
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